#!/usr/bin/R

read.params=function() {
	foo=scan(file=paste(paramDir,'pondSim.params',sep=''),what=double(),flush=TRUE,quiet=TRUE)
        delta.t<<-foo[1]
        grav<<-foo[2]
        mgd.to.cfs<<-foo[3]
	p.bay<<-foo[4]
	con.bay<<-foo[5]
	con.allen<<-foo[6]
	con.hauser<<-foo[7]
	ff.allen<<-foo[8]*mgd.to.cfs
	p.cl<<-foo[9]
        p.area<<-scan(file=paste(paramDir,'pond.area',sep=''),what=double(),flush=TRUE,quiet=TRUE)
        p.min.surf<<-scan(file=paste(paramDir,'pond.min.surf',sep=''),what=double(),flush=TRUE,quiet=TRUE)
        p.max.surf<<-scan(file=paste(paramDir,'pond.max.surf',sep=''),what=double(),flush=TRUE,quiet=TRUE)
        p.names<<-scan(file=paste(paramDir,'pond.names',sep=''),what=character(),flush=TRUE,quiet=TRUE)
        con.type<<-scan(file=paste(paramDir,'con.type',sep=''),what=integer(),flush=TRUE,quiet=TRUE)
        con.upstream<<-scan(file=paste(paramDir,'con.upstream',sep=''),what=integer(),flush=TRUE,quiet=TRUE)
        con.dnstream<<-scan(file=paste(paramDir,'con.dnstream',sep=''),what=integer(),flush=TRUE,quiet=TRUE)
        con.width<<-scan(file=paste(paramDir,'con.width',sep=''),what=double(),flush=TRUE,quiet=TRUE)
        con.l.up<<-scan(file= paste(paramDir,'con.length.up',sep=''),what=double(),flush=TRUE,quiet=TRUE)
        con.l.dn<<-scan(file= paste(paramDir,'con.length.dn',sep=''),what=double(),flush=TRUE,quiet=TRUE)
        con.d.up<<-scan(file= paste(paramDir,'con.diam.up',sep=''),what=double(),flush=TRUE,quiet=TRUE)
        con.d.dn<<-scan(file= paste(paramDir,'con.diam.dn',sep=''),what=double(),flush=TRUE,quiet=TRUE)
        con.f.up<<-scan(file= paste(paramDir,'con.friction.up',sep=''),what=double(),flush=TRUE,quiet=TRUE)
        con.f.dn<<-scan(file= paste(paramDir,'con.friction.dn',sep=''),what=double(),flush=TRUE,quiet=TRUE)
        con.names<<-scan(file=paste(paramDir,'con.names',sep=''),what=character(),flush=TRUE,quiet=TRUE)
        con.crest<<-scan(file=paste(paramDir,'con.crest',sep=''),what=double(),flush=TRUE,quiet=TRUE)
        con.top<<-scan(file=paste(paramDir,'con.top',sep=''),what=double(),flush=TRUE,quiet=TRUE)
        con.crest<<-con.top-con.crest #in the param file, values are distance from top to crest of weir
	n.ponds<<-length(p.names)
	n.cons<<-length(con.names)
	n.weirs<<-length(con.crest)
}
###############################################################################################
###############################################################################################
read.pumps = function() {
	n.pump.bins=10
	n.pumps=n.cons-n.weirs
	con.pump=array(0,dim=c(n.cons,n.pump.bins))
	#use 999 as an arbitrary large number to represent high head/flow values
	con.pump.criteria=array(999,dim=c(n.cons,n.pump.bins+1))
	for (con in 1:n.pumps) {
		b=n.weirs+con
		temp=read.table(file=paste(paramDir,'con.pump',sep=''),skip=2*con-2,nrows=1,sep='/',comment.char='#')
		con.pump[b,1:length(temp)]=as.numeric(temp) #fill in length(temp) pumping rates
		#fill in the remaining pumping rates with the highest value
		con.pump[b,(length(temp)+1):(n.pump.bins)]=as.numeric(temp)[length(temp)] 
		temp=read.table(file=paste(paramDir,'con.pump',sep=''),skip=2*con-1,nrows=1,sep='/',comment.char='#')
		#fill in pumping criteria
		con.pump.criteria[b,1:length(temp)]=as.numeric(temp)
	}	
	con.pump<<-con.pump*mgd.to.cfs #in the param file, values are given in mgd, but we need cfs
	con.pump.criteria <<- con.pump.criteria
	con.pump.criteria[con.type==4] <<- con.pump.criteria[con.type==4]*mgd.to.cfs
}
stop_ = function(){
    write(file=paste(path,'modelRunning',sep=''),'1');
    stop();
}
###############################################################################################
###############################################################################################
read.input = function() {
        p.surf<<-scan(file=paste(inputDir,'pond.surf',sep=''),
		what=double(),flush=TRUE,quiet=TRUE)
	plant.inflow<<-scan(file=paste(inputDir,'inflow',sep=''),
		what=double(),flush=TRUE,quiet=TRUE)*mgd.to.cfs #daily values (mgd, convert to cfs)
	PP<<-scan(file=paste(inputDir,'precip',sep=''),
		what=double(),flush=TRUE,quiet=TRUE)/12/24/60/60 #(file in inches/day, convert to ft/sec)
	ET<<-scan(file=paste(inputDir,'evap',sep=''),
		what=double(),flush=TRUE,quiet=TRUE)/12/24/60/60 #(file in inches/day, convert to ft/sec)
	SSE<<-scan(file=paste(inputDir,'sealevel',sep=''),
		what=double(),flush=TRUE,quiet=TRUE)
	t.horizon<<-(length(plant.inflow))*24*60*60 #t.horizon in seconds
	
	if (length(SSE)<length(plant.inflow)*4*24) {
		write(file=logfile,'***INSUFFICIENT SEALEVEL DATA, SIMULATION STOPPED***')
		stop_()
	}
	if (length(PP)<(t.horizon/60/60/24)) {
		write(file=logfile,'***INSUFFICIENT PRECIP DATA, SIMULATION STOPPED***')
		stop_()
        }
	if (length(ET)<(t.horizon/60/60/24)) {
                write(file=logfile,'***INSUFFICIENT EVAP DATA, SIMULATION STOPPED***')
		stop_()
        }
}
###############################################################################################
###############################################################################################
find.q.rectangular.weir = function(h.pond,h.friction,b) {
	my.q=0
	reduced=FALSE
	if ((h.pond)[1]>con.crest[b] ) {
		if (con.dnstream[b]==7 | h.pond[2]<h.pond[1]) { #not submerged, forward flow
			if ((h.pond-h.friction)[1]>con.crest[b]) { #(head-friction) above weir
				#francis formula:
				my.q=10/3*con.width[b]*((h.pond-h.friction)[1]-con.crest[b])^(3/2)
			}
			else if ((h.pond-h.friction)[1]<=con.crest[b]) { #(head-friction) below weir
				my.q=((h.pond[1]-con.crest[b])*(con.d.up[b]^5)*(pi^2)*grav/
				(8*con.f.up[b]*con.l.up[b]))^(1/2)	
				reduced=TRUE
        		}
		}
	}
	if (b==22) { 
		ht=1
		if (h.pond[1]>con.crest[b] & h.pond[2]>con.crest[b]) {
			if (h.pond[1]>h.pond[2] &  (h.pond[1]-h.pond[2])<(ht+1)) { 
				my.q=10/3*con.width[b]*((h.pond[1]-h.pond[2]))^(3/2) }
			if (h.pond[1]>h.pond[2] & (h.pond[1]-h.pond[2])>(ht+1)) { 
                	        my.q=10/3*con.width[b]*(ht+1)^(3/2) }
			if (h.pond[2]>h.pond[1] & (h.pond[2]-h.pond[1])<ht) { 
				my.q=(-10)/3*con.width[b]*((h.pond[2]-h.pond[1]))^(3/2) }
			if (h.pond[2]>h.pond[1] & (h.pond[2]-h.pond[1])>ht) { 
                	        my.q=(-10)/3*con.width[b]*(ht)^(3/2) }
		}
		if (h.pond[1]>con.crest[b] & h.pond[2]<=con.crest[b]) {
                        if (h.pond[1]<(con.crest[b]+ht+1)) {
                        	my.q=10/3*con.width[b]*(h.pond[1]-con.crest[b])^(3/2) }
                        if (h.pond[1]>(con.crest[b]+ht+1)) {
                        	my.q=10/3*con.width[b]*(ht+1)^(3/2) }
		}
		if (h.pond[2]>con.crest[b] & h.pond[1]<=con.crest[b]) {
                        if (h.pond[2]<(con.crest[b]+ht)) {
                                my.q=(-10)/3*con.width[b]*(h.pond[2]-con.crest[b])^(3/2) }
                        if (h.pond[2]>(con.crest[b]+ht)) {
                                my.q=(-10)/3*con.width[b]*(ht)^(3/2) }
                }
	}
        return(list(my.q,reduced))
}
###############################################################################################
###############################################################################################
find.q.vnotch.weir = function(h.pond,h.friction,b) {
	#note: the following only applies to 'right-angle' v-notch weirs
	#here, we assume no frictional losses through pipes
        my.q=0
        reduced=FALSE
        if ((h.pond)[1]>con.crest[b] & (h.pond)[2]<=con.crest[b]) { #not submerged, forward flow
		if ((h.pond-h.friction)[1]>con.crest[b]) { #(head-friction) above weir
			my.q=2.5*(h.pond[1]-h.friction[1]-con.crest[b])^(5/2)
		}
		else if ((h.pond-h.friction)[1]<con.crest[b])  { #(head-friction) below weir
			my.q=((h.pond[1]-con.crest[b])*(con.d.up[b]^5)*(pi^2)*grav/
                        (8*con.f.up[b]*con.l.up[b]))^(1/2)
                        reduced=TRUE
		}
	}
	max.q=2.5*(1/2*con.width[b])^(5/2)
	return(list(my.q,reduced))
}
####################################################################################
####################################################################################
find.flow = function(my.h) {
	#initialize local variables
	tol.flow=0.1 #max difference in flow rate between iterations
	q.old=q.new=array(0,dim=n.cons)

        for (b in 1:n.cons) {
		upstream=con.upstream[b] #index of upstream pond
		downstream=con.dnstream[b] #index of downstream pond
		up.dn = c(upstream,downstream) #indices only
                len=c(con.l.up[b],con.l.dn[b]) #pipe lengths
                friction=c(con.f.up[b],con.f.dn[b]) #pipe friction factors
                diam=c(con.d.up[b],con.d.dn[b]) #diameter of pipes
		h.local=my.h[up.dn] #heads immediately up/downstream from box
	
                
		if (con.type[b]==1) { #box is a rectangular weir
			
			#first, compute flowrate given heads upstream/downstream
			q.new[b]=find.q.rectangular.weir(h.local,c(0,0),b)[[1]]
			
			#now find head loss in pipes, re-calc flowrate, iterate until convergence
			i=0
			while(abs(abs(q.new[b])-abs(q.old[b]))>tol.flow) {
				i=i+1
				q.old[b]=q.new[b]
				h.friction = 8*friction*len*q.old[b]^2/(diam^5*pi^2*grav)
				h.friction[h.friction>abs(h.local[1]-h.local[2])]=abs(h.local[1]-h.local[2])
				foo=find.q.rectangular.weir(h.local,h.friction,b)
				q.new[b]=foo[[1]]
				if (foo[[2]]) break #friction is controlling factor, stop_ iteration
				if (i==10000) {
					write(file=logfile,append=T,paste('10000 iterations on box=',b,'time=',t))
					write(file=logfile,append=T,paste('SIMULATION STOPPED',date()))
					stop_()
				}
			}
                }
		else if (con.type[b]==2) { #box is a vnotch weir

                        #first, compute flowrate given heads upstream/downstream
                        q.new[b]=find.q.vnotch.weir(h.local,c(0,0),b)[[1]]

                        #now find head loss in pipes, re-calc flowrate, iterate until convergence
                        i=0
                        while(abs(abs(q.new[b])-abs(q.old[b]))>tol.flow) {
                                i=i+1
                                q.old[b]=q.new[b]
                                h.friction = 8*friction*len*q.old[b]^2/(diam^5*pi^2*grav)
                                h.friction[h.friction>abs(h.local[1]-h.local[2])]=abs(h.local[1]-h.local[2])
                                foo=find.q.vnotch.weir(h.local,h.friction,b)
                                q.new[b]=foo[[1]]
                                if (foo[[2]]) break
                                if (i==100) {
                                        write(file=logfile,append=T,paste('10000 iterations on box=',b,'time=',t))
                                        write(file=logfile,append=T,paste('SIMULATION STOPPED',date()))
                                        stop_()
                                }
                        }
                }
		else if (con.type[b]==3) { #box is a pump, regulated by head in upstream body
			for (crit in 1:length(con.pump[b,])) {
				#chose pump rate, based on head-criteria
				if (h.local[1]>=con.pump.criteria[b,crit] & h.local[1]<con.pump.criteria[b,crit+1]) {
					q.new[b]=con.pump[b,crit]	
				}
			}
		}
		else if (con.type[b]==4) { #box is a pump, regulated by inflow to upstream body
                        for (crit in 1:length(con.pump[b,])) {
                                #chose pump rate, based on flow-criteria
                                if (pond.inflow[con.upstream[b]]>=con.pump.criteria[b,crit] & 
					pond.inflow[con.upstream[b]]<con.pump.criteria[b,crit+1]) {
                                        q.new[b]=con.pump[b,crit]
                                }
                        }
                }	
		else if (con.type[b]==5) { #box is a pump, regulated by hours and offset
			my.begin=con.pump.criteria[b,2]*60*60
			my.end=(con.pump.criteria[b,1]+con.pump.criteria[b,2])*60*60
			if (con.pump.criteria[b,1] == -1) {
				q.new[b]=con.pump[b,1]
			}
			else if (total.t[t]>=my.begin & total.t[t]<=my.end) {
				q.new[b]=con.pump[b,1]
			}
		}	
	}
        
	#smooth the pump flow from hauser to chlorine
	tol.hauser=.05
	if (t>1) {
		if (abs(q.new[con.hauser]-ff[con.hauser,t-1]) > tol.hauser) {
			if (q.new[con.hauser]>ff[con.hauser,t-1]) { q.new[con.hauser]=ff[con.hauser,t-1]+tol.hauser }
			else { q.new[con.hauser]=ff[con.hauser,t-1]-tol.hauser }
		}
	}
		
	#if (t>20) { q.new[con.hauser]=mean(c(ff[con.hauser,(t-20):(t-1)],q.new[con.hauser])) }
	#else if (t>10) { q.new[con.hauser]=mean(c(ff[con.hauser,(t-10):(t-1)],q.new[con.hauser])) }
	#else if (t>1) { q.new[con.hauser]=mean(c(ff[con.hauser,1:(t-1)],q.new[con.hauser])) }
	
	#store flow rate in my.inflow/my.outflow vectors
	my.inflow=my.outflow=array(0,dim=n.ponds)
	for (b in 1:n.cons) {
		my.outflow[con.upstream[b]]=my.outflow[con.upstream[b]]+q.new[b]
        	my.inflow[con.dnstream[b]]=my.inflow[con.dnstream[b]]+q.new[b]
		ff[b,t]<<-q.new[b]
	}
	
	return(list(my.inflow,my.outflow))
} #end function find.flow()
#################################################################
#################################################################
cl.head.to.flow.rect=function(h) {
        flow=0
        if(h>con.crest[con.bay]) {
                flow=flow+10/3*con.width[con.bay]*(h-con.crest[con.bay])^(3/2)
        }
        return(flow)
}
cl.head.to.flow.vnotch=function(foo) {
        h=foo[[1]]; crest=foo[[2]]
        flow=0
        if(h>crest) flow=flow+2.5*(h-crest)^(5/2)
        return(flow)
}
cl.head.to.flow=function(foo) {
        h=foo[[1]]; crest=foo[[2]]
        flow=cl.head.to.flow.rect(h)+cl.head.to.flow.vnotch(list(h,crest))
        return(flow)
}
cl.flow.to.head=function(foo) {
        q=foo[[1]]; crest=foo[[2]]
        temp.func=function(h) {
                return(abs(cl.head.to.flow(list(h,crest))-q))
        }
        h=as.numeric(optimize(temp.func,lower=(crest-2),upper=(crest+2))[1])
        return(h)
}
cl.flow.to.crest=function(q) {
        temp.func=function(crest) {
                h=cl.flow.to.head(list(7,crest))
                #iterate until head and flow agree
                diff.flow=10
                to.vnotch=to.rect=0
                ii=0
                while (diff.flow>.01) {
                        ii=ii+1
                        to.vnotch = cl.head.to.flow.vnotch(list(h,crest))
                        to.rect = cl.head.to.flow.rect(h)
                        if ((to.vnotch+to.rect)>q) { h = h-.001 }
                        else { h=h+.001 }
                        diff.flow=abs(to.rect+to.vnotch-q)
                        if (ii>1000) break
                }
                #require flow to allen = flow to bay
                return(abs(to.vnotch-ff.allen))
        }
        crest=as.numeric(optimize(temp.func,lower=5,upper=10)[1])
        return(crest)
}
#################################################################
#################################################################
update.data = function() {
	total.t<<-total.t[1:(t-1)];
	ff<<-ff[,1:(t-1)]
	h<<-h[,1:(t-1)]
	t<<-t-1
	save.image(paste(path,'pondSim.RData',sep=''))
}
#################################################################
#################################################################
make.bargraph=function() {	
	h[12,h[12,]<0]=0
	h.max=array()
	for (ii in 1:n.ponds) {
		h.max[ii]=max(h[ii,])
	}
	pdf(file=paste(path,'pondSim.bargraph.pdf',sep=''),height=6,width=12)
	barplot(names.arg=p.names,height=p.max.surf,col='black',ylab='Elevation (ft above MSL)',cex.names=.9) #berms
	barplot(height=h[,1],col='gray',add=T) #initial pond heights
	barplot(height=h.max,col='darkred',add=T,width=.7,space=.686) #max pond heights reached during simulation
	barplot(height=h[,t],col='blue',add=T,width=.7,space=.686) #final pond heights
	dev.off()

    lines=readLines(paste(path,'bargraph.tex',sep=''))
    line=grep ("%Enter",lines)
    line=line+2
    c=1
    while (c<=length(h[,1])){
        currentLine=strsplit(lines[line]," ")[[1]]
        currentLine=currentLine[currentLine!=""]
        lines[line]=paste(currentLine[1],'{',h[,1][c],'}')
        line=line+1
        currentLine=strsplit(lines[line]," ")[[1]]
        currentLine=currentLine[currentLine!=""]
        lines[line]=paste(currentLine[1],'{',h.max[c],'}')
        line=line+1
        currentLine=strsplit(lines[line]," ")[[1]]
        currentLine=currentLine[currentLine!=""]
        lines[line]=paste(currentLine[1],'{',h[,t][c],'}')
        line=line+1
        c=c+1
        line=line+1
    }

    write(file=paste(path,'Hopefully.tex',sep=''),lines)
    system(paste("pdflatex -output-directory=/tmp/pondFun ",path,"Hopefully.tex",sep=''))

    write(file=paste(inputDir,"pond.surfmax",sep=''),ncolumns=1,h.max)
    write(file=paste(inputDir,"pond.surffinal",sep=''),ncolumns=1,h[,t])
    write(file=paste(inputDir,"pond.surfinitial",sep=''),ncolumns=1,h[,1])
    
}

################################################################
################################################################
